/* copyright (c) 1996, 1997, 1998, 1999 William R. Pearson and the
   U. of Virginia */

/* $Name: fasta34t11d4 $ - $Id: comp_mlib.c,v 1.8 2002/06/24 13:04:17 wrp Exp $ */

/*
 * Concurrent read version
 *
 *	Feb 20, 1998 modifications for prss3
 *
 *	December, 1998 - DNA searches are now down with forward and reverse
 *			 strands
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <time.h>

#include <limits.h>
#include <float.h>
#include <math.h>

#ifdef UNIX
#include <unistd.h>
#include <sys/types.h>
#include <signal.h>
#endif

#ifdef __MWERKS__
#include <Types.h>
#include <StandardFile.h>
StandardFileReply freply;
Point wpos;
int tval;
char prompt[256];
#include <sioux.h>
#endif

#include "defs.h"
#include "mm_file.h"

#include "mw.h"			/* defines beststr */
#include "structs.h"		/* mngmsg, libstruct */
#include "param.h"		/* pstruct, thr_str, buf_head, rstruct */

#define XTERNAL
#include "uascii.h"

char *mp_verstr="34t11";

/********************************/
/* global variable declarations */
/********************************/
char gstring2[MAX_STR];                  /* string for label */
char gstring3[MAX_STR];
char hstring1[MAX_STR];

extern int max_workers;

extern int nsfnum;
extern int sfnum[];
#ifdef SUPERFAMNUM
extern int sfn_cmp(int *q, int *s);
extern int nsfnum_n;
extern int sfnum_n[];
#endif

#ifdef __MWERKS__
FSSpec ouSpec, q0Spec, q1Spec;
OSErr error;
extern int PgmDID;
int ccommand(char ***);	/* MSL */
void InitEvent();	/* checkevent.c */
void ChkEvent();	/* checkevent.c */
#endif

/********************************/
/* extern variable declarations */
/********************************/
extern char *prog_func;		/* function label */
extern char *verstr, *iprompt0, *iprompt1, *iprompt2, *refstr;

/********************************/
/*extern function declarations  */
/********************************/

/* open sequence file (getseq.c) */
extern int getseq(char *filen, unsigned char *seq, int maxs, char *libstr,
		   int *seq_type, int seq_ext, long *sq0ff);

#ifndef PRSS
struct lmf_str *openlib(char *, int, int *, int);
#else
struct lmf_str *openlib(char *, int, int *, struct mngmsg);
#endif

void closelib(struct lmf_str *m_fptr);

void irand(int);

/* check for DNA sequence (nmgetlib.c) */
extern int scanseq(unsigned char *seq, int n, char *str);
extern void revcomp(unsigned char *seq, int n, int *c_nt);

/* initialize environment (doinit.c) */
extern void initenv (int argc, char **argv, struct mngmsg *m_msg,
		     struct pstruct *ppst, unsigned char **aa0);

/* print timing information */
extern void ptime (FILE *, time_t);

/* read a sequence (nmgetlib.c) */
extern int (*getlib) (unsigned char *seq, int maxs,
		      char *libstr, 
		      int n_libstr,
		      FSEEK_T *libpos,
		      int *lcont,
		      long *l_off);

#define QGETLIB (q_file_p->getlib)
#define GETLIB (m_file_p->getlib)

/* calculation functions */
extern void
init_work(unsigned char *aa0, int n0, struct pstruct *ppst, void **f_arg);
#ifndef COMP_THR
extern void
do_work(unsigned char *aa0, int n0, unsigned char *aa1, int n1, int frame, 
	struct pstruct *ppst, void *f_str, struct rstruct *rst);
#endif
extern void
close_work(unsigned char *aa0, int n0, struct pstruct *ppst, void **f_arg);
extern void
get_param (struct pstruct *pstr, char *pstring1, char *pstring2);

#ifdef COMP_THR
#ifndef PRSS
void
save_best(struct buf_head *cur_buf, struct pstruct pst, int nitt1,
	  FILE *fdata, int *, struct hist_str *, void **);
#else
void
save_best(struct buf_head *cur_buf, struct pstruct pst, int nitt1,
	  FILE *fdata, int *, struct hist_str *, void **, int *, int *);
#endif
#endif

/* statistics functions */
extern void
process_hist(struct stat_str *sptr, int nstat, struct pstruct pst,
	     struct hist_str *hist, void **);
extern void addhistz(float, struct hist_str *); /* scaleswn.c */
void selectbestz(struct beststr **, int, int );
extern double zs_to_E(float, int, int, long, struct db_str);
extern float (*find_zp)(int score, int length, double comp, void *);

#ifndef COMP_THR
extern int shuffle(unsigned char *, unsigned char *, int);
extern int wshuffle(unsigned char *, unsigned char *, int, int, int *);
#endif

extern void irand(int);
extern void set_db_size(int, struct db_str *, struct hist_str *);

/* display functions */
#ifndef PRSS
extern void
showbest (FILE *fp, unsigned char **aa0, unsigned char *aa1,
	  int maxn, struct beststr **bestp_arr, int nbest,
	  int qlib, struct mngmsg *m_msg,struct pstruct pst,
	  struct db_str db, void **f_str, char *gstring2);

#else
extern void
showbest (FILE *fp, int *s_score, int maxn, int qlib,
	  struct mngmsg *m_msg,struct pstruct pst,
	  struct db_str db, void *, char *gstring2);
#endif

extern void
showalign (FILE *fp, unsigned char **aa0, unsigned char *aa1,
	   int maxn, struct beststr **bestp_arr, int nbest,
	   int qlib, struct mngmsg m_msg,struct pstruct pst,
	   void **f_str, char *gstring2);

/* misc functions */
void h_init(struct pstruct *);		/* doinit.c */
void last_init(struct mngmsg *, struct pstruct *); /* initfa/sw.c */
void s_abort(char *, char *);		/* compacc.c */

/* initfa/sw.c */
void resetp(struct mngmsg *, struct pstruct *); 

void libchoice(char *lname, int, struct mngmsg *); /* lib_sel.c */
void libselect(char *lname, struct mngmsg *);	/* lib_sel.c */
void query_parm(struct mngmsg *, struct pstruct *); /* initfa/sw.c */
void selectbestz(struct beststr **, int, int);
void sortbestz(struct beststr **, int);
void sortbeste(struct beststr **, int);
void sortbest(struct beststr **, int, int);
void prhist(FILE *, struct mngmsg, struct pstruct, 
	    struct hist_str hist, int nstats, struct db_str, char *);
void printsum(FILE *, struct db_str db);

FILE *outfd;			/* Output file */

/* this information is global for fsigint() */
extern time_t s_time ();		/* fetches time */
time_t tstart, tscan, tprev, tdone, total_scan, total_disp;	/* Timing */
time_t tdstart, tddone, time();

static struct db_str qtt = {0l, 0l, 0};

#ifdef COMP_THR
/***************************************/
/* thread global variable declarations */
/***************************************/

/* functions for getting/sending buffers to threads (thr_sub.c) */
extern void init_thr(int , struct thr_str **, struct mngmsg, struct pstruct *,
		     unsigned char *, int);
extern void start_thr(void);
extern void get_rbuf(struct buf_head **cur_buf, int max_wor_buf);
extern void put_rbuf(struct buf_head *cur_buf, int max_work_buf);
extern void put_rbuf_done(int nthreads, struct buf_head *cur_buf, 
			  int max_work_buf);
#undef XTERNAL
#include "thr.h"
struct buf_head buf_list[NUM_WORK_BUF];
#endif

/* these variables must be global for comp_thr.c so that savebest()
   can use them */

static struct beststr 
    *best,		/* array of best scores */
    *bestp,
    **bestp_arr;	/* array of pointers */
static int nbest;	/* number of best scores */

static struct stat_str *stats;	/* array of scores for statistics */
static int nstats;
static float zbestcut;	/* cut off for best z-score */
static int bestfull;			/* index for selectbest() */
static int dohist=0;			/* flag for z-value processing */

main (int argc, char *argv[])
{
  unsigned char *aa0[6], *aa1, *aa1ptr, *aa1s;
  int n1;

  int *n1tot_ptr=NULL, *n1tot_cur;
  int n1tot_cnt=0;
  int n1tot_v;

  long l_off;
  char lib_label[MAX_FN];
  char q_bline[MAX_STR];
  long qseek;
  int qlib;
  struct lmf_str *m_file_p, *q_file_p;

  int t_best, t_rbest;		/* best score of two/six frames */
  int i_score;
#ifdef PRSS
  int s_score[3];
  int s_n1;
#endif

  struct pstruct pst;
  void *f_str[6];		/* different f_str[]'s for different
				   translation frames, or forward,reverse */
  int have_f_str=0;

#ifdef COMP_THR
  long ntbuff;
  int max_buf_cnt, ave_seq_len, buf_siz;
  int max_work_buf;
  struct buf_head *cur_buf;
  struct buf_str *cur_buf_p;
  int nseq;
  struct thr_str *work_info[MAX_WORKERS];
#endif

  struct mngmsg m_msg;		/* Message from host to manager */
  int iln, itt;			/* index into library names */
  char rline[MAX_FN];
  int t_quiet;

  struct a_struct *aln_d_base, *aln_d_ptr;

  struct rstruct  rst;		/* results structure */
  struct rstruct  rrst;		/* results structure for shuffle*/
  int i,j;

  FILE *fdata=NULL;		/* file for full results */
  char libstr[21];		/* string for labeling full results */
  
  int leng;			/* leng is length of the descriptive line */
  int maxn = MAXTOT;		/* max size of the query sequence */
  FSEEK_T lmark;		/* seek into library of current sequence */
  int lcont, ocont, maxt;	/* continued sequence */
  int igncnt=0;			/* count for ignoring sequences warning */
  int ieven=0;			/* tmp for wshuffle */
  float zscore;			/* tmp value */
  char *bp;			/* general purpose string ptr */
  void fsigint();
  
  /* Initialization */

#ifdef __MWERKS__
  SIOUXSettings.showstatusline=FALSE;
  SIOUXSettings.asktosaveonclose=TRUE;
  SIOUXSettings.autocloseonquit=TRUE;
  
  argc = ccommand(&argv);

  InitEvent();
  /* GetVol((unsigned char *)prompt,&ouvRef); */
  error=HGetVol(NULL,&ouSpec.vRefNum, &ouSpec.parID);
  if (error != noErr) {
  	fprintf(stderr," cannot get current directory\n");
  	exit(1);
  	}
  wpos.h=50; wpos.v=100;
#endif

  m_msg.quiet = 1;

  /* first initialization routine - nothing is known */
  h_init(&pst);
  
  f_str[0] = f_str[1] = NULL;

  aa0[0] = NULL;
  /* second initialiation - get commmand line arguments */
  initenv (argc, argv, &m_msg, &pst,&aa0[0]);

#ifdef COMP_THR
  /* now have max_workers - allocate work_info[] */
  for (i=0; i<max_workers; i++) {
    if ((work_info[i]=
	 (struct thr_str *)calloc(sizeof(1),sizeof(struct thr_str)))==NULL) {
      fprintf(stderr, " cannot allocate work_info[%d]\n",i);
      exit(1);
    }
  }
#endif
#ifndef PRSS
  /* label library size limits */
  if (m_msg.n1_low > 0 && m_msg.n1_high < BIGNUM) 
    sprintf(lib_label,"library (range: %d-%d)",m_msg.n1_low,m_msg.n1_high);
  else if (m_msg.n1_low > 0) 
    sprintf(lib_label,"library (range: >%d)",m_msg.n1_low);
  else if (m_msg.n1_high < BIGNUM)
    sprintf(lib_label,"library (range: <%d)",m_msg.n1_high);
  else
    strncpy(lib_label,"library",sizeof(lib_label));
#else
  sprintf(lib_label,"shuffled sequence");
#endif
  lib_label[sizeof(lib_label)-1]='\0';

#ifdef __MWERKS__
  if (m_msg.quiet) { SIOUXSettings.asktosaveonclose=FALSE; }
#endif

  total_scan = total_disp = 0;
  tstart = tscan = s_time();
  tdstart = time(NULL);

  /* Allocate space for the query and library sequences */
  if (aa0[0]==NULL) {
    if ((aa0[0] = (unsigned char *)malloc((MAXTOT+1)*sizeof(unsigned char)))
	== NULL)
      s_abort ("Unable to allocate query sequence", "");
    *aa0[0]=0;
    aa0[0]++;
  }
  aa0[5]=aa0[4]=aa0[3]=aa0[2]=aa0[1]=aa0[0];

  /* make room for random sequence */
  if ((aa1s = (unsigned char *)malloc((MAXTOT+1) * sizeof (char))) == NULL) {
      s_abort ("Unable to allocate shuffled library sequence", "");
      *aa1s=0;
      aa1s++;
  }
  irand(0);

  if (m_msg.markx == 6) {
    fprintf(stdout,"<html>\n<head>\n<title>%s Results\n",
	    prog_func);
    fprintf(stdout,"</title>\n</head>\n<body>\n<pre>\n");
  }

  fputs(iprompt0,stdout);
  fprintf(stdout," %s%s\n",verstr,refstr);

  /* Query library */
  if (m_msg.tname[0] == '\0') {
      if (m_msg.quiet == 1)
	s_abort("Query sequence undefined","");
#ifndef __MWERKS__           
    l1:	fputs (iprompt1, stdout);
      fflush  (stdout);
      if (fgets (m_msg.tname, MAX_FN, stdin) == NULL)
	s_abort ("Unable to read query library name","");
      if ((bp=strchr(m_msg.tname,'\n'))!=NULL) *bp='\0';
      if (m_msg.tname[0] == '\0') goto l1;
#else
    NIntroDlog(PgmDID,verstr,"\0","\0","\0");	

  l1:	SFileDlog(iprompt1,&freply);
    if (freply.sfGood==TRUE) {
      PtoCstr((unsigned char *)freply.sfFile.name);
      strcpy(m_msg.tname,(char *)freply.sfFile.name);
	  q1Spec.vRefNum = q0Spec.vRefNum = freply.sfFile.vRefNum;
	  q1Spec.parID = q0Spec.parID = freply.sfFile.parID;
	  HSetVol(NULL,q0Spec.vRefNum,q0Spec.parID);
    }
    else exit(0);
#endif
  }
  

  /* Open query library */
  if ((q_file_p=
#ifndef PRSS
       openlib(m_msg.tname, m_msg.qdnaseq,qascii,!m_msg.quiet))==NULL)
#else
       openlib(m_msg.tname, m_msg.qdnaseq,qascii,m_msg))==NULL)
#endif
    {
      s_abort(" cannot open library ",m_msg.tname);
    }
  qlib = 0;

  resetp (&m_msg, &pst);

  /* get library file names */

#ifdef __MWERKS__
    HSetVol(NULL,ouSpec.vRefNum,ouSpec.parID);
#endif

#ifndef PRSS
    if (strlen (m_msg.lname) == 0) {
      if (m_msg.quiet == 1) s_abort("Library name undefined","");
      libchoice(m_msg.lname,sizeof(m_msg.lname),&m_msg);
    }
  
    libselect(m_msg.lname, &m_msg);
#else
    if (strlen (m_msg.lname) == 0) {
      if (m_msg.quiet == 1) s_abort("Shuffle sequence undefined","");
    l2:    fputs(iprompt2,stdout);
      fflush(stdout);
      if (fgets (m_msg.lname, MAX_FN, stdin) == NULL)
	s_abort ("Unable to read shuffle file name","");
      if ((bp=strchr(m_msg.lname,'\n'))!=NULL) *bp='\0';
      if (m_msg.lname[0] == '\0') goto l2;
    }
    m_msg.lbnames[0]= m_msg.lname;
    m_msg.nln = 1;
    m_msg.nshow = 0;
#endif

    /* Get additional parameters here */
    if (!m_msg.quiet) query_parm (&m_msg, &pst);
  
    last_init(&m_msg, &pst);

  /* Allocate space for saved scores */
    if ((best = 
	 (struct beststr *)malloc((MAXBEST+1)*sizeof(struct beststr)))==NULL)
      s_abort ("Cannot allocate best struct","");
    if ((bestp_arr = 
	 (struct beststr **)malloc((MAXBEST+1)*sizeof(struct beststr *)))==NULL)
      s_abort ("Cannot allocate bestp_arr","");
  
  /* Initialize bestp_arr */
    for (nbest = 0; nbest < MAXBEST+1; nbest++)
      bestp_arr[nbest] = &best[nbest];
    best++; bestp_arr++;
    best[-1].score[0]=best[-1].score[1]=best[-1].score[2]= INT_MAX;
    best[-1].zscore=FLT_MAX;
  
    if ((stats =
	 (struct stat_str *)malloc((MAXSTATS)*sizeof(struct stat_str)))==NULL)
      s_abort ("Cannot allocate stats struct","");

    outfd = stdout;

#ifdef UNIX
    /* set up signals now that input is done */
    signal(SIGHUP,SIG_IGN);
#endif

#ifdef COMP_THR
    /* Set up buffers for reading the library:

     We will start by using a 2 Mbyte buffer for each worker.  For
     proteins, that means 5,000 sequences of length 400 (average).
     For DNA, that means 2,000 sequences of length 1000.  At the
     moment, those are good averages.
  */

    if (m_msg.ldnaseq==1) {
      max_buf_cnt = MAX_NT_BUF;
      ave_seq_len = AVE_NT_LEN;
    }
    else {
      max_buf_cnt = MAX_AA_BUF;
      ave_seq_len = AVE_AA_LEN;
    }

    /* however - buffer sizes should be a function of the number of
     workers so that all the workers are kept busy.  Assuming a 10,000
     entry library is the smallest we want to schedule, then
  */

    if (max_buf_cnt > 10000/max_workers) 
      max_buf_cnt = 10000/(2*max_workers);

    max_buf_cnt /= m_msg.thr_fact;

    /* finally, max_work_buf should be mod 6 for tfasta */
    max_buf_cnt -= (max_buf_cnt % 6);

    max_work_buf = 2*max_workers;

  /* allocate space for library buffers and results */

    buf_siz=max_buf_cnt*ave_seq_len;
    if (buf_siz < 3*maxn) buf_siz = 3*maxn;
    for (i=0; i<max_work_buf; i++) {
      if ((buf_list[i].buf =(struct buf_str *)calloc((size_t)(max_buf_cnt+1),
						     sizeof(struct buf_str)))
	  ==NULL) {
	fprintf(stderr," cannot allocate buffer struct %d %d\n",i,max_buf_cnt+1);
	exit(1);
      }
      buf_list[i].buf_cnt=0;
      buf_list[i].have_results=0;
      if ((buf_list[i].start =
	   (unsigned char *)calloc((size_t)(buf_siz),sizeof(unsigned char)))
	  ==NULL) {
	fprintf(stderr," cannot allocate buffer %d\n",i);
	exit(1);
      }
      reader_buf[i] = &buf_list[i];
    }

    /* initialization of global variables for threads/buffers */

    num_worker_bufs = 0;
    num_reader_bufs = max_work_buf;
    reader_done = 0;
    worker_buf_workp = 0;
    worker_buf_readp = 0;
    reader_buf_workp = 0;
    reader_buf_readp = 0;

    start_thread = 1;	/* keeps threads from starting */
#endif

  printf("Query library %s vs %s library\n", m_msg.tname,m_msg.lname);
  printf("searching %s library\n",m_msg.lbnames[0]);

  m_msg.pstat_void = NULL;

  qtt.length = 0l;
  qtt.entries = qtt.carry = 0;

  /* Fetch first sequence */
  lcont = 0;
  while (1) {

    nbest = 0;
    zbestcut = -FLT_MAX;
    nstats = dohist= 0;

    pst.n0 = m_msg.n0 = 
      QGETLIB (aa0[0], MAXTST, m_msg.qtitle, sizeof(m_msg.qtitle), &qseek, &lcont,q_file_p,&m_msg.sq0off);

    if (m_msg.n0 <= 0) break;

    m_msg.db.length = 0l;
    m_msg.db.entries = m_msg.db.carry = 0;

    strncpy(m_msg.libstr,m_msg.qtitle,sizeof(m_msg.libstr));
    m_msg.libstr[sizeof(m_msg.libstr)-1]='\0';
    if ((bp=strchr(m_msg.libstr,' '))!=NULL) *bp='\0';

    qlib++;
    qtt.length += m_msg.n0;

    if (qlib <= 0) {
      fprintf(stderr," no sequences found in query library\n");
      exit(1);
    }

#ifdef SUPERFAMNUM
    m_msg.nqsfnum = nsfnum;
    for (i=0; i <= nsfnum & i<10; i++) m_msg.qsfnum[i] = sfnum[i];
    m_msg.nqsfnum_n = nsfnum_n;
    for (i=0; i <= nsfnum_n & i<10; i++) m_msg.qsfnum_n[i] = sfnum_n[i];
#endif

    maxn = MAXTOT - m_msg.n0 -2;	/* maxn = max library sequence space */
    /* reduce maxn if requested */
    if (m_msg.maxn > 0 && m_msg.maxn < maxn) maxn = m_msg.maxn;

#ifndef COMP_THR
    /* serial does revcomp() here, threaded does it in work_thr.c */
    if (m_msg.qframe == 2) {
      aa0[1] = aa0[0] + m_msg.n0 + 2;
      *aa0[1] = '\0';
      aa0[1]++;
      maxn -= m_msg.n0+2;
      memcpy(aa0[1],aa0[0],m_msg.n0);
      revcomp(aa0[1],m_msg.n0,&pst.c_nt[0]);
      /*
	for (i=0; i<m_msg.n0; i++) {
	fputc(nt[aa0[1][i]],stderr);
	if (i%60==59) fputc('\n',stderr);
	}
	fputc('\n',stderr);
      */
    }
    /* set aa1 for serial - threaded points aa1 to buffer */
    aa1 = aa0[1] + m_msg.n0+2;
#else
    init_thr(max_workers, work_info, m_msg, &pst, aa0[0], max_work_buf);
#endif

    if (m_msg.qdnaseq==m_msg.ldnaseq || m_msg.qdnaseq==1) {  /* !TFAST */
      if (m_msg.n0>(MAXTST+MAXLIB)/3) {
	fprintf(stderr,
		" query sequence is too long %d %s\n",m_msg.n0,m_msg.sqnam);
	exit(1);
      }
      m_msg.loff = m_msg.n0;
      m_msg.maxt3 = maxn-m_msg.loff;
    }
    else {	/* TFAST */
      if (m_msg.n0 > MAXTST) {
	fprintf(stderr,
		" query sequence is too long %d %s\n",m_msg.n0,m_msg.sqnam);
	exit(1);
      }

      if (m_msg.n0*6 > maxn ) {	/* n0 for the complement + n0*3 for
                                   the three frames */
	if (m_msg.n0*7+2 < MAXTST+MAXLIB) { /* m_msg.n0*6 + m_msg.n0 */
	  fprintf(stderr,
		  " query sequence too long for library segment: %d - resetting to %d\n",
		  maxn,m_msg.n0*6);
	  maxn = m_msg.maxn = m_msg.n0*6;
	}
	else {
	  fprintf(stderr,
		  " query sequence too long for translated search: %d %s\n",
		  m_msg.n0,m_msg.sqnam);
	  exit(1);
	}
      }

      /* set up some constants for overlaps */

      m_msg.loff = 3*m_msg.n0;
      m_msg.maxt3 = maxn-m_msg.loff-3;
      m_msg.maxt3 -= m_msg.maxt3%3;
      m_msg.maxt3++;

      maxn = maxn - 3; maxn -= maxn%3; maxn++;
    }

    pst.maxlen = maxn;
  
    leng = (int)strlen (m_msg.libstr);
    sprintf (&m_msg.libstr[leng], " %d %s", m_msg.n0, m_msg.sqnam);
    m_msg.seqnm = 0;
  
    /* Label the output */
    if ((bp = (char *) strchr (m_msg.lname, ' ')) != NULL) *bp = '\0';
  
    if (m_msg.ltitle[0] == '\0') {
      strncpy(m_msg.ltitle,m_msg.lname,sizeof(m_msg.ltitle));
      m_msg.ltitle[sizeof(m_msg.ltitle)-1]='\0';
    }

#ifndef PRSS
    /* rline[] is a tmp string */
    if (m_msg.qdnaseq==1) {
      strncpy(rline,(m_msg.qframe==1)? " (forward-only)" : "\0",sizeof(rline));
      rline[sizeof(rline)-1]='\0';
    }
    else rline[0]='\0';

    printf("%3d>>>%.50s: %d %s%s\n", qlib,
	   m_msg.qtitle, m_msg.n0, m_msg.sqnam,
	   (m_msg.revcomp ? " (reverse complement)" : rline));
#endif

    if (outfd != stdout)
      fprintf (outfd, "Query sequence: %.50s (%d %s)%s\n vs  %.60s %s\n",
	       m_msg.tname, m_msg.n0, m_msg.sqnam, (m_msg.revcomp)?
	       " (reverse-complement)":rline,m_msg.lname,
	       lib_label);
  
    fflush (outfd);

    tprev = s_time();
  
    if (m_msg.dfile[0] && (fdata=fopen(m_msg.dfile,"w"))!=NULL)
      fprintf(fdata,"%3d\t%-50s\n",m_msg.n0,m_msg.qtitle);

    qtt.length += m_msg.n0;
    qtt.entries++;

#ifdef COMP_THR
    start_thr();

  /* now open the library and start reading */
  /* get a buffer and fill it up */
    get_rbuf(&cur_buf,max_work_buf);

    cur_buf->buf_cnt = 0;
    cur_buf->have_results = 0;
    cur_buf->buf[0].aa1b = cur_buf->start;
    ntbuff = 0;
    nseq = 0;
#else
    /* initialize the comparison function, returning f_str */
    init_work (aa0[0], m_msg.n0, &pst, &f_str[0]);
    have_f_str = 1;

    f_str[5] = f_str[4] = f_str[3] = f_str[2] = f_str[1] = f_str[0];
    if (m_msg.qframe == 2) {
      init_work ( aa0[1], m_msg.n0, &pst, &f_str[1]);
    }
#endif

    /* open the library - start the search */

    for (iln = 0; iln < m_msg.nln; iln++) {
      if ((m_msg.lb_mfd[iln] = m_file_p=
#ifndef PRSS
	   openlib(m_msg.lbnames[iln], m_msg.ldnaseq, lascii, !m_msg.quiet))
#else
	   openlib(m_msg.lbnames[iln], m_msg.ldnaseq, lascii, m_msg))
#endif
	   ==NULL) {
           fprintf(stderr," cannot open library %s\n",m_msg.lbnames[iln]);
           continue;
      }

    lcont = 0;
    ocont = 0;
    n1tot_v = n1tot_cnt = 0;
    n1tot_cur = n1tot_ptr = NULL;

    /* get next buffer to read into */
    maxt = maxn;

#ifndef COMP_THR
    aa1ptr = aa1;
#else
    /* read sequence directly into buffer */
    aa1ptr = aa1 = cur_buf->buf[nseq].aa1b;
#endif

    while ((n1=GETLIB(aa1ptr,maxt,libstr,20,&lmark,&lcont,m_file_p,&l_off))>=0) {
      if ((bp=strchr(libstr,' '))!=NULL) *bp='\0';

      /* check for a continued sequence and provide a pointer to 
	 the n1_tot array if lcont || ocont */
      n1tot_v += n1;
      if (lcont && !ocont) {	/* get a new pointer */
	if (n1tot_cnt <= 0) {
	  if ((n1tot_ptr=calloc(1000,sizeof(int)))==NULL) {
	    fprintf(stderr," cannot allocate n1tot_ptr\n");
	    exit(1);
	  }
	  else {n1tot_cnt=1000;}
	}
	n1tot_cnt--;
	n1tot_cur = n1tot_ptr++;
      }

      if (n1tot_v < m_msg.n1_low || n1tot_v > m_msg.n1_high) {
	goto loop2;
      }

      m_msg.db.entries++;
      m_msg.db.length += n1;
      if (m_msg.db.length > LONG_MAX) {
	m_msg.db.length -= LONG_MAX; m_msg.db.carry++;
      }

#ifdef DEBUG
      /* This finds most reasons for core dumps */
      if (pst.debug_lib)
	for (i=0; i<n1; i++)
	  if (aa1[i]>=pst.nsqx) 
	    {fprintf(stderr,
		     "%s residue[%d/%d] %d range (%d) lcont/ocont: %d/%d\n%s\n",
		     libstr,i,n1,aa1[i],pst.nsqx,lcont,ocont,aa1ptr+i);
	    aa1[i]=0;
	    n1=i-1;
	    break;
	    }
#endif

      /* don't count long sequences more than once */
      if (aa1!=aa1ptr) {n1 += m_msg.loff; m_msg.db.entries--;}

#ifdef PROGRESS
#ifdef UNIX
      if (!m_msg.quiet) 
#endif
	if (m_msg.db.entries % 200 == 199) {
	  fputc('.',stderr);
	  if (m_msg.db.entries % 10000 == 9999) fputc('\n',stderr);
	  else if (m_msg.db.entries % 1000 == 999) fputc(' ',stderr);
#ifdef __MWERKS__
	  fflush(stderr);	/* this is a bug */
#endif
	}
#endif
#ifdef __MWERKS__
      ChkEvent();
#endif

      if (n1<=1) {
	/*	if (igncnt++ <10)
		fprintf(stderr,"Ignoring: %s\n",libstr);
	*/
	goto loop2;
      }

      /* if COMP_THR - fill and empty buffers */
#ifdef COMP_THR
      ntbuff += n1+1;

      for (itt=m_msg.revcomp; itt<=m_msg.nitt1; itt++) {

	cur_buf->buf_cnt++;
	cur_buf_p = &(cur_buf->buf[nseq++]);
	cur_buf_p->n1  = n1;
	cur_buf_p->n1tot_p = n1tot_cur;
	cur_buf_p->lseek = lmark;
	cur_buf_p->cont = ocont+1;
	cur_buf_p->m_file_p = (void *)m_file_p;
	cur_buf_p->frame = itt;
	strncpy(cur_buf_p->libstr,libstr,12);
	cur_buf_p->libstr[12]='\0';
	cur_buf_p->nsfnum = nsfnum;
#ifdef SUPERFAMNUM
	if ((cur_buf_p->sfnum[0]=sfnum[0])>0 &&
	    (cur_buf_p->sfnum[1]=sfnum[1])>0 &&
	    (cur_buf_p->sfnum[2]=sfnum[2])>0 &&
	    (cur_buf_p->sfnum[3]=sfnum[3])>0 &&
	    (cur_buf_p->sfnum[4]=sfnum[4])>0 &&
	    (cur_buf_p->sfnum[5]=sfnum[5])>0 &&
	    (cur_buf_p->sfnum[6]=sfnum[6])>0 &&
	    (cur_buf_p->sfnum[7]=sfnum[7])>0 &&
	    (cur_buf_p->sfnum[8]=sfnum[8])>0 &&
	    (cur_buf_p->sfnum[9]=sfnum[9])>0) ;
#endif

	/* this assumes that max_buf_cnt is guaranteed %6=0 so that
	   additional pointers to the same buffer can be used 
	   nseq now points to next buffer
	*/

	cur_buf->buf[nseq].aa1b = cur_buf->buf[nseq-1].aa1b;
      } /* for (itt .. */

      /* make a copy of the overlap (threaded only) */
      if (lcont) {
	memcpy(aa1s,&aa1[n1-m_msg.loff],m_msg.loff);
      }

      /* if the buffer is filled */
      if (nseq >= max_buf_cnt || ntbuff >= buf_siz - maxn) {
	/* provide filled buffer to workers */
	put_rbuf(cur_buf,max_work_buf);

	/* get an empty buffer to fill */
	get_rbuf(&cur_buf,max_work_buf);

	/* "empty" buffers have results that must be processed */
	if (cur_buf->buf_cnt && cur_buf->have_results) {
#ifndef PRSS
	  save_best(cur_buf,pst,m_msg.nitt1,fdata,m_msg.qsfnum,&m_msg.hist,
		    &m_msg.pstat_void);
#else
	  save_best(cur_buf,pst,m_msg.nitt1,fdata,m_msg.qsfnum,&m_msg.hist,
		    &m_msg.pstat_void,s_score,&s_n1);
#endif
	}

	/* now the buffer is truly empty, fill it up */
	cur_buf->buf_cnt = 0;
	cur_buf->have_results = 0;
	/* point the first aa1 ptr to the buffer start */
	aa1=cur_buf->buf[0].aa1b = cur_buf->start;
	ntbuff = 0;
	nseq=0;
      }
      else {	/* room left in current buffer, increment ptrs */
	aa1=cur_buf->buf[nseq].aa1b = cur_buf->buf[nseq-1].aa1b+n1+1;
      }
#else /* if !COMP_THR - do a bunch of searches */

      /* t_best and t_rbest are used to save the best score or shuffled
	 score from all the frames */

      t_best = t_rbest = -1;
      for (itt=m_msg.revcomp; itt<=m_msg.nitt1; itt++) {

	rst.score[0] = rst.score[1] = rst.score[2] = 0;
	do_work (aa0[itt], m_msg.n0, aa1, n1, itt, &pst, f_str[itt], &rst);

	if (rst.score[pst.score_ix] > t_best) {
	  t_best = rst.score[pst.score_ix];
	}

	if (fdata) {
	  fprintf(fdata,"%-12s %5d %6d %d %.5lf %.5lf %4d %4d %4d %8ld\n",
		  libstr,sfn_cmp(m_msg.qsfnum,sfnum),n1,itt,
		  rst.comp,rst.H,
		  rst.score[0],rst.score[1],rst.score[2],
		  lmark);
	  fflush(fdata);
	}

#ifdef PRSS
	if (lmark==0) {
	  s_score[0] = rst.score[0];
	  s_score[1] = rst.score[1];
	  s_score[2] = rst.score[2];
	  s_n1 = n1;
	  m_msg.db.entries=0;
	  m_msg.db.length=0;
	  continue;
	}
	t_best = t_rbest = rst.score[pst.score_ix];
#else
	if (pst.zsflag >= 10) {
	  if (pst.zs_win > 0) wshuffle(aa1,aa1s,n1,pst.zs_win,&ieven);
	  else shuffle(aa1,aa1s,n1);
	  do_work (aa0[itt], m_msg.n0, aa1s, n1, itt, &pst, f_str[itt], &rrst);
	  if (rrst.score[pst.score_ix] > t_rbest) {
	    t_rbest = rrst.score[pst.score_ix];
	  }
	}
#endif
	i_score = rst.score[pst.score_ix];
	if (nstats < MAXSTATS) {
	  if (itt==m_msg.nitt1) {	/* all frames are done */
	    stats[nstats].n1 = n1;	/* save the best score */
	    stats[nstats].comp =  rst.comp;
	    stats[nstats].H = rst.H;
	    if (pst.zsflag >=10) t_best = t_rbest;
	    stats[nstats++].score = t_best;
	    t_best = t_rbest = -1;	/* reset t_rbest, t_best */
	  }
	}
	else if (!dohist && pst.zsflag >= 0) {
  /* this process_hist is called in the middle of a search */
	  process_hist(stats,nstats,pst,&m_msg.hist,&m_msg.pstat_void);
	  for (i=0; i<nbest; i++) {
	    bestp_arr[i]->zscore = 
	      (*find_zp)(bestp_arr[i]->score[pst.score_ix],
			 bestp_arr[i]->n1, bestp_arr[i]->comp, m_msg.pstat_void);
	  }
	  dohist = 1;
	}

	if (dohist) { /* dohist > 0 => nstats >= MAXSTATS */
	  zscore=(*find_zp)(i_score,n1, rst.comp,m_msg.pstat_void);

	  if (itt == m_msg.nitt1) {
	    if (pst.zsflag >= 10) t_best = t_rbest;
	    
	    addhistz((*find_zp)(t_best,n1, rst.comp, m_msg.pstat_void),
		     &m_msg.hist);
	    t_best = t_rbest = -1;
	  }
	}
	else zscore = (float) i_score;

	if (zscore > zbestcut ) {
	  if (nbest >= MAXBEST) {
	    bestfull = nbest-MAXBEST/4;
	    selectbestz(bestp_arr,bestfull-1,nbest);
	    zbestcut = bestp_arr[bestfull-1]->zscore;
	    nbest = bestfull;
	  }

	  bestp = bestp_arr[nbest++];
	  bestp->score[0] = rst.score[0];
	  bestp->score[1] = rst.score[1];
	  bestp->score[2] = rst.score[2];
	  bestp->comp =  rst.comp;
	  bestp->H = rst.H;
	  bestp->zscore = zscore;
	  bestp->lseek = lmark;
	  memcpy(bestp->libstr,libstr,MAX_UID);
	  bestp->cont = ocont+1;
	  bestp->m_file_p = m_file_p;
	  bestp->n1 = n1;
	  bestp->n1tot_p=n1tot_cur;
	  bestp->frame = itt;
	  bestp->nsfnum = nsfnum;
#ifdef SUPERFAMNUM
	  if ((bestp->sfnum[0]=sfnum[0])>0 &&
	      (bestp->sfnum[1]=sfnum[1])>0 &&
	      (bestp->sfnum[2]=sfnum[2])>0 &&
	      (bestp->sfnum[3]=sfnum[3])>0 &&
	      (bestp->sfnum[4]=sfnum[4])>0 &&
	      (bestp->sfnum[5]=sfnum[5])>0 &&
	      (bestp->sfnum[6]=sfnum[6])>0 &&
	      (bestp->sfnum[7]=sfnum[7])>0 &&
	      (bestp->sfnum[8]=sfnum[8])>0 &&
	      (bestp->sfnum[9]=sfnum[9])>0) ;
#endif
	}
      }
#endif

    loop2: 
      if (lcont) {
	maxt = m_msg.maxt3;
#ifndef COMP_THR
	memcpy(aa1,&aa1[n1-m_msg.loff],m_msg.loff);
#else
	memcpy(aa1,aa1s,m_msg.loff);
#endif
	aa1ptr= &aa1[m_msg.loff];
	ocont = lcont;
      }
      else {
	maxt = maxn;
	aa1ptr=aa1;
	if (ocont) *n1tot_cur = n1tot_v;
	ocont = 0;
	n1tot_v = 0;
	n1tot_cur = NULL;
      }
    } /* end while((n1=getlib())) */
  } /* end iln=1..nln */

  /* all done */

#ifdef COMP_THR
  /* check last buffers for any results */
  put_rbuf_done(max_workers,cur_buf,max_work_buf);

  for (i=0; i < num_reader_bufs; i++) {
    reader_buf_readp = (reader_buf_readp+1)%(max_work_buf);
    if (reader_buf[reader_buf_readp]->buf_cnt > 0 && 
	reader_buf[reader_buf_readp]->have_results) {
#ifndef PRSS
      save_best( reader_buf[reader_buf_readp],pst,m_msg.nitt1,fdata,m_msg.qsfnum,
		 &m_msg.hist, &m_msg.pstat_void);
#else
      save_best( reader_buf[reader_buf_readp],pst,m_msg.nitt1,fdata,m_msg.qsfnum,
		 &m_msg.hist, &m_msg.pstat_void,s_score, &s_n1);
#endif
    }
  }
#endif

#ifdef PROGRESS
#ifdef UNIX
  if (!m_msg.quiet) {
#endif
    if (m_msg.db.entries >= 200) fprintf(stderr," Done!\n");
    /* this is a bug */
#ifdef __MWERKS__
    fflush(stderr);
#endif
  }
#endif

  m_msg.nbr_seq = m_msg.db.entries;
  get_param(&pst, gstring2,gstring3);

  /* *************************** */
  /* analyze the last results    */
  /* *************************** */
    
#ifndef PRSS
  if (!dohist) {
    if ((nstats < 20 || pst.zsflag < 0) && pst.zsflag != 3) {
      pst.zsflag = -1;
    }
    else {

  /* this process_hist is called when everything is done, and does not
     require additional addhistz()'s */

      process_hist(stats,nstats,pst,&m_msg.hist, &m_msg.pstat_void);
      for (i = 0; i < nbest; i++) {
	bestp_arr[i]->zscore = (*find_zp)(bestp_arr[i]->score[pst.score_ix],
					  bestp_arr[i]->n1, bestp_arr[i]->comp,
					  m_msg.pstat_void);
      }
    }
  }
#else
  /* this process_hist is called when everything is done, only for PRSS */
  if (pst.zsflag < 10) pst.zsflag += 10;
  process_hist(stats,nstats,pst,&m_msg.hist, &m_msg.pstat_void);
  for (i = 0; i < nbest; i++) {
    bestp_arr[i]->zscore = (*find_zp)(bestp_arr[i]->score[pst.score_ix],
				      bestp_arr[i]->n1,bestp_arr[i]->comp,
				      m_msg.pstat_void);
  }
#endif

  if (pst.zsflag >=0 ) {
    sortbestz (bestp_arr, nbest);
    if (pst.zdb_size <= 1) pst.zdb_size = m_msg.db.entries;
    for (i=0; i<nbest; i++) {
      bestp_arr[i]->escore=zs_to_E(bestp_arr[i]->zscore,bestp_arr[i]->n1,pst.dnaseq,
				   pst.zdb_size,m_msg.db);
    }
    sortbeste(bestp_arr, nbest);
  }
  else sortbest(bestp_arr,nbest,pst.score_ix);

  m_msg.db.entries = m_msg.db.entries;
  m_msg.db.length = m_msg.db.length;
  m_msg.db.carry = m_msg.db.carry;

#ifdef PRSS
  gettitle(m_msg.lname,m_msg.ltitle,sizeof(m_msg.ltitle));
  printf("%.50s - %s %d %s%s\n vs %.60s - %s shuffled sequence\n",
	 m_msg.tname, m_msg.qtitle,m_msg.n0, m_msg.sqnam,
	 (m_msg.revcomp ? " (reverse complement)" : "\0"),
	 m_msg.lname,m_msg.ltitle);
#endif

  prhist (stdout, m_msg, pst, m_msg.hist, nstats, m_msg.db, gstring2);

  tscan = s_time();
  printf (" Scan time: ");
  total_scan += tscan-tprev;
  ptime(stdout,tscan-tprev);
  printf ("\n");

 l3:
  if (!m_msg.quiet) {
    printf("Enter filename for results [%s]: ", m_msg.outfile);
    fflush(stdout);
  }

  rline[0]='\0';
  if (!m_msg.quiet && fgets(rline,sizeof(rline),stdin)==NULL) goto end_l;
  if ((bp=strchr(rline,'\n'))!=NULL) *bp = '\0';
  if (rline[0]!='\0') strncpy(m_msg.outfile,rline,sizeof(m_msg.outfile));
  if (m_msg.outfile[0]!='\0') {
    if (outfd==NULL && (outfd=fopen(m_msg.outfile,"w"))==NULL) {
      fprintf(stderr," could not open %s\n",m_msg.outfile);
      if (!m_msg.quiet) goto l3;
      else goto l4;
    }
#ifdef __MWERKS__
    HSetVol(NULL,q0Spec.vRefNum,q0Spec.parID);
#endif
    /* label the output */
    if (m_msg.markx == 6) {
      fprintf(outfd,"<html>\n<head>\n<title>%s Results\n",
	      prog_func);
      fprintf(outfd,"</title>\n</head>\n<body>\n<pre>\n");
    }

    fputs(iprompt0,outfd);
    fprintf(outfd," %s%s\n",verstr,refstr);

    fprintf(outfd," %s%s, %d %s\n vs %s %s\n",
	    m_msg.tname, (m_msg.revcomp ? "-" : "\0"), m_msg.n0,
	    m_msg.sqnam, m_msg.ltitle, lib_label);

    prhist(outfd,m_msg,pst,m_msg.hist, nstats, m_msg.db, gstring2);
  }

 l4:   
  if (m_msg.markx == 6) {
    fprintf(outfd,"</pre>\n<p>\n<hr>\n<p>\n");
  }

#ifdef PRSS
  showbest (stdout, s_score, s_n1, qtt.entries, &m_msg, pst, m_msg.db, 
	    m_msg.pstat_void,gstring2);
  if (outfd != stdout) {
    showbest (outfd, s_score, maxn, qtt.entries, &m_msg, pst, m_msg.db, 
	      m_msg.pstat_void,gstring2);
  }
#else
  /* code from p2_complib.c to pre-calculate -m 9 alignment info -
       requires -q with -m 9 */

  if (m_msg.quiet || m_msg.markx==9) {

    /* to determine how many sequences to re-align (either for
       do_opt() or calc_id() we need to modify m_msg.mshow to get
       the correct number of alignments */

    if (m_msg.mshow_flg != 1 && pst.zsflag >= 0) {
      for (i=0; i<nbest && bestp_arr[i]->escore< m_msg.e_cut; i++) {}
      m_msg.mshow = i;
    }

    if (m_msg.mshow <= 0) { /* no results to display */
      fprintf(outfd,"!! No sequences with E() < %f\n",m_msg.e_cut);
      m_msg.nshow = 0;
      goto end_l;
    }
  }

#ifdef COMP_THR
  /* before I call showbest/showalign, I need init_work() to get an f_str */
  init_work(aa0[0],m_msg.n0,&pst,&f_str[0]);
  have_f_str = 1;
  f_str[5] = f_str[4] = f_str[3] = f_str[2] = f_str[1] =  f_str[0];

  if (m_msg.qframe == 2) {
    if ((aa0[1]=(unsigned char *)calloc((size_t)m_msg.n0+2,
					sizeof(unsigned char)))==NULL) {
      fprintf(stderr," cannot allocate aa0[1][%d] for alignments\n",
	      m_msg.n0+2);
    }
    *aa0[1]='\0';
    aa0[1]++;
    memcpy(aa0[1],aa0[0],m_msg.n0+1);
    revcomp(aa0[1],m_msg.n0,&pst.c_nt[0]);
    init_work(aa0[1],m_msg.n0,&pst,&f_str[1]);
  }

  /* I also need a "real" aa1 */
  aa1 = buf_list[0].start;
#endif

  showbest (stdout, aa0, aa1, maxn, bestp_arr, nbest, qtt.entries, &m_msg, pst,
	    m_msg.db, f_str, gstring2);

  if (outfd != stdout) {
    t_quiet = m_msg.quiet;
    m_msg.quiet = -1;	/* should guarantee 1..nbest shown */
    showbest (outfd, aa0, aa1, maxn, bestp_arr, nbest, qtt.entries, &m_msg, pst,
	      m_msg.db, f_str, gstring2);
    m_msg.quiet = t_quiet;
  }
#endif

  if (m_msg.nshow > 0) {
    rline[0]='N';
    if (!m_msg.quiet){
      printf(" Display alignments also? (y/n) [n] "); fflush(stdout);
      if (fgets(rline,sizeof(rline),stdin)==NULL) goto end_l;
    }
    else rline[0]='Y';

    if (toupper((int)rline[0])=='Y') {
      if (!m_msg.quiet) {
	printf(" number of alignments [%d]? ",m_msg.nshow);
	fflush(stdout);
	if (fgets(rline,sizeof(rline),stdin)==NULL) goto end_l;
	if (rline[0]!=0) sscanf(rline,"%d",&m_msg.nshow);
	m_msg.ashow=m_msg.nshow;
      }

      if (m_msg.markx == 5 || m_msg.markx==6 || m_msg.markx==9) {
	fprintf(outfd,"\n>>>%s%s, %d %s vs %s library\n",
		m_msg.tname,(m_msg.revcomp ? "_rev":"\0"), m_msg.n0,
		m_msg.sqnam,m_msg.lname);
      }

      if (m_msg.markx==10) {
	fprintf(outfd,"\n>>>%s%s, %d %s vs %s library\n",
		m_msg.tname,(m_msg.revcomp ? "-":"\0"), m_msg.n0, m_msg.sqnam,
		m_msg.lname);
	fprintf(outfd,"; mp_name: %s\n",argv[0]);
	fprintf(outfd,"; mp_ver: %s\n",mp_verstr);
	fprintf(outfd,"; mp_argv:");
	for (i=0; i<argc; i++)
	  fprintf(outfd," %s",argv[i]);
	fputc('\n',outfd);
	fputs(gstring3,outfd);
	fputs(hstring1,outfd);
      }

      showalign (outfd, aa0, aa1, maxn, bestp_arr, nbest, qtt.entries, m_msg, pst,
		 f_str, gstring2);

      fflush(outfd);
    }
  }

 end_l:
#ifdef COMP_THR
  for (i=0; i<max_work_buf; i++) {
      buf_list[i].buf_cnt=0;
      buf_list[i].have_results=0;
  }

  num_worker_bufs = 0;
  num_reader_bufs = max_work_buf;
  reader_done = 0;
  worker_buf_workp = 0;
  worker_buf_readp = 0;
  reader_buf_workp = 0;
  reader_buf_readp = 0;

  start_thread = 1;	/* stop thread from starting again */
#endif

  if (have_f_str) {
    if (f_str[1]!=f_str[0]) {
      close_work (aa0[1], m_msg.n0, &pst, &f_str[1]);
    }
    close_work (aa0[0], m_msg.n0, &pst, &f_str[0]);
    have_f_str = 0;
  }

  for (iln=0; iln < m_msg.nln; iln++) {
    if (m_msg.lb_mfd[iln]!=NULL) closelib(m_msg.lb_mfd[iln]);
  }

  tddone = time(NULL);
  tdone = s_time();
  total_disp += tdone-tscan;
  fflush(outfd);

    if (fdata) {
      fprintf(fdata,"/** %s **/\n",gstring2);
      fprintf(fdata,"%3ld%-50s\n",qtt.entries-1,m_msg.libstr);
      fflush(fdata);
    }
  }

  if (m_msg.markx==10)
      fprintf(outfd,">>><<<\n");

    tdone = s_time();
#ifndef PRSS
    if ( m_msg.markx == 6) fprintf(outfd,"<p><pre>\n"); 
    printsum(outfd, m_msg.db);
    if ( m_msg.markx == 6) fprintf(outfd,"</pre>\n");
#endif
    if (m_msg.markx == 6) fprintf(outfd,"</body></html>\n");
    if (outfd!=stdout) printsum(stdout,m_msg.db);

    exit(0);
}   /* End of main program */

void
printsum(FILE *fd, struct db_str ntt)
{
  double db_tt;
  char tstr1[26], tstr2[26];

  strncpy(tstr1,ctime(&tdstart),sizeof(tstr1));
  strncpy(tstr2,ctime(&tddone),sizeof(tstr1));
  tstr1[24]=tstr2[24]='\0';

  /* Print timing to output file as well */
  fprintf(fd, "\n\n%ld residues in %ld query   sequences\n", qtt.length, qtt.entries);
  if (ntt.carry == 0) 
    fprintf(fd, "%ld residues in %ld library sequences\n", ntt.length, ntt.entries);
  else {
    db_tt = (double)ntt.carry*(double)LONG_MAX + (double)ntt.length;
    fprintf(fd, "%.0f residues in %ld library sequences\n", db_tt, ntt.entries);
  }

#ifndef COMP_THR
  fprintf(fd," Scomplib [%s]\n start: %s done: %s\n",mp_verstr,tstr1,tstr2);
#else
  fprintf(fd," Tcomplib [%s] (%d proc)\n start: %s done: %s\n", mp_verstr,
    max_workers,tstr1,tstr2);
#endif
  fprintf(fd," Scan time: ");
  ptime(fd, total_scan);
  fprintf (fd," Display time: ");
  ptime (fd, total_disp);
  fprintf (fd,"\n");
  fprintf (fd, "\nFunction used was %s [%s]\n", prog_func,verstr);
}

void fsigint()
{
  struct db_str db;

  db.entries = db.length = db.carry = 0;
  tdone = s_time();
  tddone = time(NULL);

  printf(" /*** interrupted ***/\n");
  if (outfd!=stdout) fprintf(outfd,"/*** interrupted ***/\n");
  fprintf(stderr,"/*** interrupted ***/\n");

  printsum(stdout,db);
  if (outfd!=stdout) printsum(outfd,db);

  exit(1);
}

#ifdef COMP_THR
#ifndef PRSS
void save_best(struct buf_head *cur_buf, struct pstruct pst, int nitt1,
	       FILE *fdata, int *qsfnum, struct hist_str *histp,
	       void **pstat_voidp)
#else
void save_best(struct buf_head *cur_buf, struct pstruct pst, int nitt1,
	       FILE *fdata, int *qsfnum, struct hist_str *histp,
	       void **pstat_voidp, int *s_score, int *s_n1)
#endif
{
  float zscore;
  register int i_score;
  struct buf_str *p_rbuf, *cur_buf_p;
  int i, t_best, t_rbest, tm_best, t_n1, sc_ix;

  sc_ix = pst.score_ix;

  cur_buf_p = cur_buf->buf;
  
  t_best = t_rbest = -1;

  while (cur_buf->buf_cnt--) { /* count down the number of results */
    p_rbuf = cur_buf_p++;	/* step through the results buffer */

    i_score = p_rbuf->rst.score[sc_ix];
    /* need to look for frame 0 if TFASTA, then save stats at frame 6 */
    if (fdata) {
      fprintf(fdata,"%-12s %5d %6d %d %.5lf %.5lf %4d %4d %4d %8ld\n",
	      p_rbuf->libstr,
#ifdef SUPERFAMNUM
	      sfn_cmp(qsfnum,p_rbuf->sfnum),
#else
	      0,
#endif
	      p_rbuf->n1,p_rbuf->frame,p_rbuf->rst.comp,p_rbuf->rst.H,
	      p_rbuf->rst.score[0],p_rbuf->rst.score[1],p_rbuf->rst.score[2],
	      p_rbuf->lseek);
    }

#ifdef PRSS
    if (p_rbuf->lseek==0) {
      s_score[0] = p_rbuf->rst.score[0];
      s_score[1] = p_rbuf->rst.score[1];
      s_score[2] = p_rbuf->rst.score[2];
      *s_n1 = p_rbuf->n1;
      continue;
    }
#endif

    if (i_score > t_best) {
      tm_best = t_best = i_score;
      t_n1 = p_rbuf->n1;
    }
    if (pst.zsflag >= 10 && p_rbuf->r_score > t_rbest)
      t_rbest = p_rbuf->r_score;

    /* statistics done for best score of set */

    if (nstats < MAXSTATS ) {
      if (p_rbuf->frame == nitt1) {
	stats[nstats].n1 = t_n1;
	stats[nstats].comp = p_rbuf->rst.comp;
	stats[nstats].H = p_rbuf->rst.H;
	if (pst.zsflag >= 10) {
	  tm_best = t_rbest;
	  t_rbest = -1;
	}
	stats[nstats++].score = tm_best;
	t_best = -1;
      }
    }
    else if (!dohist && pst.zsflag > 0) {
      process_hist(stats,nstats,pst, histp, pstat_voidp);
      dohist = 1;
      for (i=0; i<nbest; i++)
	  bestp_arr[i]->zscore = 
	    (*find_zp)(bestp_arr[i]->score[pst.score_ix],bestp_arr[i]->n1,bestp_arr[i]->comp,
		       *pstat_voidp);
    }

    /* best saved for every score */
    if (dohist) {

      zscore=(*find_zp)(i_score, p_rbuf->n1,(float)p_rbuf->rst.comp,
			*pstat_voidp);

      if (p_rbuf->frame == nitt1) {
	addhistz((*find_zp)(t_best, p_rbuf->n1, (float) p_rbuf->rst.comp,
			    *pstat_voidp), histp);
	t_best = t_rbest = -1;
      }
    }
    else zscore = (float) i_score;

    if (zscore > zbestcut) {
      if (nbest >= MAXBEST) {
	bestfull = nbest-MAXBEST/4;
	selectbestz(bestp_arr,bestfull-1,nbest);
	zbestcut = bestp_arr[bestfull-1]->zscore;
	nbest = bestfull;
      }
      bestp = bestp_arr[nbest++];
      bestp->score[0] = p_rbuf->rst.score[0];
      bestp->score[1] = p_rbuf->rst.score[1];
      bestp->score[2] = p_rbuf->rst.score[2];
      bestp->comp = (float) p_rbuf->rst.comp;
      bestp->H = (float) p_rbuf->rst.H;
      bestp->zscore = zscore;
      bestp->lseek = p_rbuf->lseek;
      memcpy(bestp->libstr,p_rbuf->libstr,MAX_UID);
      bestp->cont = p_rbuf->cont; /* not cont+1 because incremented already */
      bestp->m_file_p = p_rbuf->m_file_p;
      bestp->n1 = p_rbuf->n1;
      bestp->n1tot_p = p_rbuf->n1tot_p;
      bestp->frame = p_rbuf->frame;
      bestp->nsfnum = p_rbuf->nsfnum;
#ifdef SUPERFAMNUM
      if ((bestp->sfnum[0] = p_rbuf->sfnum[0])>0 &&
	  (bestp->sfnum[1] = p_rbuf->sfnum[1])>0 &&
	  (bestp->sfnum[2] = p_rbuf->sfnum[2])>0 &&
	  (bestp->sfnum[3] = p_rbuf->sfnum[3])>0 &&
	  (bestp->sfnum[4] = p_rbuf->sfnum[4])>0 &&
	  (bestp->sfnum[5] = p_rbuf->sfnum[5])>0 &&
	  (bestp->sfnum[6] = p_rbuf->sfnum[6])>0 &&
	  (bestp->sfnum[7] = p_rbuf->sfnum[7])>0 &&
	  (bestp->sfnum[8] = p_rbuf->sfnum[8])>0 &&
	  (bestp->sfnum[9] = p_rbuf->sfnum[9])>0) ;
#endif
    }
  }
}
#endif
